{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "<center>\n",
    "<img src=\"../../img/ods_stickers.jpg\">\n",
    "## Открытый курс по машинному обучению\n",
    "<center>Автор материала: Юрий Кашницкий\n",
    "    \n",
    "Материал распространяется на условиях лицензии [Creative Commons CC BY-NC-SA 4.0](https://creativecommons.org/licenses/by-nc-sa/4.0/). Можно использовать в любых целях (редактировать, поправлять и брать за основу), кроме коммерческих, но с обязательным упоминанием автора материала."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <center>Тема 4. Линейные модели классификации и регрессии\n",
    "## <center>Часть 2. Логистическая регрессия и метод максимального правдоподобия "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Линейный классификатор"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Основная идея линейного классификатора заключается в том, что признаковое пространство может быть разделено гиперплоскостью на две полуплоскости, в каждой из которых прогнозируется одно из двух значений целевого класса. \n",
    "Если это можно сделать без ошибок, то обучающая выборка называется *линейно разделимой*.\n",
    "\n",
    "<img src=\"../../img/logit.png\">\n",
    "\n",
    "Мы уже знакомы с линейной регрессией и методом наименьших квадратов. Рассмотрим задачу бинарной классификации, причем метки целевого класса обозначим \"+1\" (положительные примеры) и \"-1\" (отрицательные примеры).\n",
    "Один из самых простых линейных классификаторов получается на основе регрессии вот таким образом:\n",
    "\n",
    "$$\\Large a(\\textbf{x}) = \\text{sign}(\\textbf{w}^{\\text{T}}\\textbf x),$$\n",
    "\n",
    "где\n",
    " - $\\textbf{x}$ – вектор признаков примера (вместе с единицей);\n",
    " - $\\textbf{w}$ – веса в линейной модели (вместе со смещением $w_0$);\n",
    " - $\\text{sign}(\\bullet)$ – функция \"сигнум\", возвращающая знак своего аргумента;\n",
    " - $a(\\textbf{x})$ – ответ классификатора на примере $\\textbf{x}$.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Логистическая регрессия как линейный классификатор"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Логистическая регрессия является частным случаем линейного классификатора, но она обладает хорошим \"умением\" – прогнозировать вероятность $p_+$ отнесения примера $\\textbf{x}_\\text{i}$ к классу \"+\":\n",
    "$$\\Large p_+ = \\text P\\left(y_i = 1 \\mid \\textbf{x}_\\text{i}, \\textbf{w}\\right) $$\n",
    "\n",
    "Прогнозирование не просто ответа (\"+1\" или \"-1\"), а именно *вероятности* отнесения к классу \"+1\" во многих задачах является очень важным бизнес-требованием. Например, в задаче кредитного скоринга, где традиционно применяется логистическая регрессия, часто прогнозируют вероятность невозврата кредита ($p_+$). Клиентов, обратившихся за кредитом, сортируют по этой предсказанной вероятности (по убыванию), и получается скоркарта — по сути, рейтинг клиентов от плохих к хорошим. Ниже приведен игрушечный пример такой скоркарты. \n",
    "    <img src='../../img/toy_scorecard.png' width=60%>\n",
    "\n",
    "Банк выбирает для себя порог $p_*$ предсказанной вероятности невозврата кредита (на картинке – $0.15$) и начиная с этого значения уже не выдает кредит. Более того, можно умножить предсказнную вероятность на выданную сумму и получить матожидание потерь с клиента, что тоже будет хорошей бизнес-метрикой."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Итак, мы хотим прогнозировать вероятность $p_+ \\in [0,1]$, а пока умеем строить линейный прогноз с помощью МНК: $b(\\textbf{x}) = \\textbf{w}^\\text{T} \\textbf{x} \\in \\mathbb{R}$. Каким образом преобразовать полученное значение в вероятность, пределы которой – [0, 1]? Очевидно, для этого нужна некоторая функция $f: \\mathbb{R} \\rightarrow [0,1].$ В модели логистической регрессии для этого берется конкретная функция: $\\sigma(z) = \\frac{1}{1 + \\exp^{-z}}$. И сейчас разберемся, каковы для этого предпосылки. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# отключим всякие предупреждения Anaconda\n",
    "import warnings\n",
    "\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sigma(z):\n",
    "    return 1.0 / (1 + np.exp(-z))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xx = np.linspace(-10, 10, 1000)\n",
    "plt.plot(xx, [sigma(x) for x in xx])\n",
    "plt.xlabel(\"z\")\n",
    "plt.ylabel(\"sigmoid(z)\")\n",
    "plt.title(\"Sigmoid function\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Обозначим $P(X)$ вероятностью происходящего события $X$. Тогда отношение вероятностей $OR(X)$ определяется из $\\frac{P(X)}{1-P(X)}$, а это — отношение вероятностей того, произойдет ли событие или не произойдет. Очевидно, что вероятность и отношение шансов содержат одинаковую информацию. Но в то время как $P(X)$ находится в пределах от 0 до 1, $OR(X)$ находится в пределах от 0 до $\\infty$.\n",
    "\n",
    "Если вычислить логарифм $OR(X)$ (то есть называется логарифм шансов, или логарифм отношения вероятностей), то легко заметить, что $\\log{OR(X)} \\in \\mathbb{R}$. Его то мы и будем прогнозировать с помощью МНК.\n",
    "\n",
    "Посмотрим, как логистическая регрессия будет делать прогноз $p_+ = \\text{P}\\left(y_i = 1 \\mid \\textbf{x}_\\text{i}, \\textbf{w}\\right)$ (пока считаем, что веса $\\textbf{w}$ мы как-то получили (т.е. обучили модель), далее разберемся, как именно). \n",
    "\n",
    "**Шаг 1.** Вычислить значение $w_{0}+w_{1}x_1 + w_{2}x_2 + ... = \\textbf{w}^\\text{T}\\textbf{x}$. (уравнение $\\textbf{w}^\\text{T}\\textbf{x} = 0$ задает гиперплоскость, разделяющую примеры на 2 класса);\n",
    "\n",
    "\n",
    "**Шаг 2.** Вычислить логарифм отношения шансов: $ \\log(OR_{+}) =  \\textbf{w}^\\text{T}\\textbf{x}$.\n",
    "\n",
    "**Шаг 3.** Имея прогноз шансов на отнесение к классу \"+\" – $OR_{+}$, вычислить $p_{+}$ с помощью простой зависимости:\n",
    "\n",
    "$$\\Large p_{+} = \\frac{OR_{+}}{1 + OR_{+}} = \\frac{\\exp^{\\textbf{w}^\\text{T}\\textbf{x}}}{1 + \\exp^{\\textbf{w}^\\text{T}\\textbf{x}}} =  \\frac{1}{1 + \\exp^{-\\textbf{w}^\\text{T}\\textbf{x}}} = \\sigma(\\textbf{w}^\\text{T}\\textbf{x})$$\n",
    "\n",
    "\n",
    "В правой части мы получили как раз сигмоид-функцию.\n",
    "\n",
    "Итак, логистическая регрессия прогнозирует вероятность отнесения примера к классу \"+\" (при условии, что мы знаем его признаки и веса модели) как сигмоид-преобразование линейной комбинации вектора весов модели и вектора признаков примера:\n",
    "\n",
    "$$\\Large p_+(x_i) = \\text{P}\\left(y_i = 1 \\mid \\textbf{x}_\\text{i}, \\textbf{w}\\right) = \\sigma(\\textbf{w}^\\text{T}\\textbf{x}_\\text{i}). $$\n",
    "\n",
    "Следующий вопрос: как модель обучается. Тут мы опять обращаемся к принципу максимального правдоподобия."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Принцип максимального правдоподобия и логистическая регрессия\n",
    "Теперь посмотрим, как из принципа максимального правдоподобия получается оптимизационная задача, которую решает логистическая регрессия, а именно, – минимизация *логистической* функции потерь. \n",
    "Только что мы увидели, что логистическая регрессия моделирует вероятность отнесения примера к классу \"+\" как \n",
    "\n",
    "$$\\Large p_+(\\textbf{x}_\\text{i}) = \\text{P}\\left(y_i = 1 \\mid \\textbf{x}_\\text{i}, \\textbf{w}\\right) = \\sigma(\\textbf{w}^\\text{T}\\textbf{x}_\\text{i})$$\n",
    "\n",
    "Тогда для класса \"-\" аналогичная вероятность:\n",
    "$$\\Large p_-(\\textbf{x}_\\text{i})  = \\text{P}\\left(y_i = -1 \\mid \\textbf{x}_\\text{i}, \\textbf{w}\\right)  = 1 - \\sigma(\\textbf{w}^\\text{T}\\textbf{x}_\\text{i}) = \\sigma(-\\textbf{w}^\\text{T}\\textbf{x}_\\text{i}) $$\n",
    "\n",
    "Оба этих выражения можно ловко объединить в одно (следите за моими руками – не обманывают ли вас):\n",
    "\n",
    "$$\\Large \\text{P}\\left(y = y_i \\mid \\textbf{x}_\\text{i}, \\textbf{w}\\right) = \\sigma(y_i\\textbf{w}^\\text{T}\\textbf{x}_\\text{i})$$\n",
    "\n",
    "Выражение $M(\\textbf{x}_\\text{i}) = y_i\\textbf{w}^\\text{T}\\textbf{x}_\\text{i}$ называется *отступом* (*margin*) классификации на объекте $\\textbf{x}_\\text{i}$ (не путать с зазором (тоже margin), про который чаще всего говорят в контексте SVM). Если он неотрицателен, модель не ошибается на объекте $\\textbf{x}_\\text{i}$, если же отрицателен – значит, класс для $\\textbf{x}_\\text{i}$  спрогнозирован неправильно. \n",
    "Заметим, что отступ определен для объектов именно обучающей выборки, для которых известны реальные метки целевого класса $y_i$.\n",
    "\n",
    "Чтобы понять, почему это мы сделали такие выводы, обратимся к геометрической интерпретации линейного классификатора. Подробно про это можно почитать в материалах Евгения Соколова – [тут](https://github.com/esokolov/ml-course-msu/blob/master/ML16/lecture-notes/Sem09_linear.pdf). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Рекомендую решить почти классическую задачу из начального курса линейной алгебры: найти расстояние от точки с радиус-вектором $\\textbf{x}_A$ до плоскости, которая задается уравнением $\\textbf{w}^\\text{T}\\textbf{x} = 0.$\n",
    "\n",
    "\n",
    "Ответ: \n",
    "$\\Large \\rho(\\textbf{x}_A, \\textbf{w}^\\text{T}\\textbf{x} = 0) = \\frac{\\textbf{w}^\\text{T}\\textbf{x}_A}{||\\textbf{w}||}$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src = '../../img/simple_linal_task.png' width=60%>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Когда получим (или посмотрим) ответ, то поймем, что чем больше по модулю выражение $\\textbf{w}^{\\text{T}}\\textbf{x}_\\text{i}$, тем дальше точка $\\textbf{x}_\\text{i}$ находится от плоскости $\\textbf{w}^{\\text{T}}\\textbf{x} = 0.$\n",
    "\n",
    "Значит, выражение $M(\\textbf{x}_\\text{i}) = y_i\\textbf{w}^{\\text{T}}\\textbf{x}_\\text{i}$ – это своего рода \"уверенность\" модели в классификации объекта $\\textbf{x}_\\text{i}$: \n",
    "\n",
    "- если отступ большой (по модулю) и положительный, это значит, что метка класса поставлена правильно, а объект находится далеко от разделяющей гиперплоскости (такой объект классифицируется уверенно). На рисунке – $x_3$.\n",
    "- если отступ большой (по модулю) и отрицательный, значит метка класса поставлена неправильно, а объект находится далеко от разделюящей гиперплоскости (скорее всего такой объект – аномалия, например, его метка в обучающей выборке поставлена неправильно). На рисунке – $x_1$.\n",
    "- если отступ малый (по модулю), то объект находится близко к разделюящей гиперплоскости, а  знак отступа определяет, правильно ли объект классифицирован.  На рисунке – $x_2$ и $x_4$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src = '../../img/margin.png' width=60%>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Теперь распишем правдоподобие выборки, а именно, вероятность наблюдать данный вектор $\\textbf{y}$ у выборки $\\textbf X$. Делаем сильное предположение: объекты приходят независимо, из одного распределения (*i.i.d.*). Тогда\n",
    "\n",
    "$$\\Large \\text{P}\\left(\\textbf{y} \\mid \\textbf X, \\textbf{w}\\right) = \\prod_{i=1}^{\\ell} \\text{P}\\left(y = y_i \\mid \\textbf{x}_\\text{i}, \\textbf{w}\\right),$$\n",
    "\n",
    "где $\\ell$ – длина выборки $\\textbf X$ (число строк).\n",
    "\n",
    "Как водится, возьмем логарифм данного выражения (сумму оптимизировать намного проще, чем произведение):\n",
    "\n",
    "$$\\Large  \\log \\text{P}\\left(\\textbf{y} \\mid \\textbf X, \\textbf{w}\\right) = \\log \\sum_{i=1}^{\\ell} \\text{P}\\left(y = y_i \\mid \\textbf{x}_\\text{i}, \\textbf{w}\\right) = \\log \\prod_{i=1}^{\\ell} \\sigma(y_i\\textbf{w}^{\\text{T}}\\textbf{x}_\\text{i})   = $$\n",
    "\n",
    "$$\\Large  = \\sum_{i=1}^{\\ell} \\log \\sigma(y_i\\textbf{w}^{\\text{T}}\\textbf{x}_\\text{i}) = \\sum_{i=1}^{\\ell} \\log \\frac{1}{1 + \\exp^{-y_i\\textbf{w}^{\\text{T}}\\textbf{x}_\\text{i}}} = - \\sum_{i=1}^{\\ell} \\log (1 + \\exp^{-y_i\\textbf{w}^{\\text{T}}\\textbf{x}_\\text{i}})$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "То есть в данном случае принцип максимизации правдоподобия приводит к минимизации выражения \n",
    "\n",
    "$$\\Large \\mathcal{L_{log}} (\\textbf X, \\textbf{y}, \\textbf{w}) = \\sum_{i=1}^{\\ell} \\log (1 + \\exp^{-y_i\\textbf{w}^{\\text{T}}\\textbf{x}_\\text{i}}).$$\n",
    "\n",
    "Это *логистическая* функция потерь, просуммированная по всем объектам обучающей выборки.\n",
    "\n",
    "Посмотрим на новую фунцию как на функцию от отступа: $L(M) = \\log (1 + \\exp^{-M})$. Нарисуем ее график, а также график 1/0 функциий потерь (*zero-one loss*), которая просто штрафует модель на 1 за ошибку на каждом объекте (отступ отрицательный): $L_{1/0}(M) = [M < 0]$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src = '../../img/logloss_margin.png' width=60%>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Картинка отражает общую идею, что в задаче классификации, не умея напрямую минимизировать число ошибок (по крайней мере, градиентными методами это не сделать – производная 1/0 функциий потерь в нуле обращается в бесконечность), мы минимизируем некоторую ее верхнюю оценку. В данном случае это логистическая функция потерь (где логарифм двоичный, но это не принципиально), и справедливо \n",
    "\n",
    "$$\\Large \\mathcal{L_{\\text{1/0}}} (\\textbf X, \\textbf{y}, \\textbf{w}) = \\sum_{i=1}^{\\ell} [M(\\textbf{x}_\\text{i}) < 0] \\leq \\sum_{i=1}^{\\ell} \\log (1 + \\exp^{-y_i\\textbf{w}^{\\text{T}}\\textbf{x}_\\text{i}}) = \\mathcal{L_{\\log}} (\\textbf X, \\textbf{y}, \\textbf{w}), $$\n",
    "\n",
    "где $\\mathcal{L_{\\text{1/0}}} (\\textbf X, \\textbf{y}, \\textbf{w})$ – попросту число ошибок логистической регрессии с весами $\\textbf{w}$ на выборке $(\\textbf X, \\textbf{y})$.\n",
    "\n",
    "То есть уменьшая верхнюю оценку $\\mathcal{L_{\\log}}$ на число ошибок классификации, мы таким образом надеемся уменьшить и само число ошибок."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### L2-регуляризация логистической функции потерь"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$L2-регуляризация$ логистической регрессии устроена почти так же, как и в случае с гребневой (Ridge регрессией). Вместо функционала $\\mathcal{L_{\\log}} (X, \\textbf{y}, \\textbf{w})$ минимизируется следующий:\n",
    "\n",
    "$$\\Large J(\\textbf X, \\textbf{y}, \\textbf{w}) = \\mathcal{L_{\\log}} (\\textbf X, \\textbf{y}, \\textbf{w}) + \\lambda |\\textbf{w}|^2$$\n",
    "\n",
    "В случае логистической регрессии принято введение обратного коэффициента регуляризации $C = \\frac{1}{\\lambda}$. И тогда решением задачи будет\n",
    "\n",
    "$$\\Large \\widehat{\\textbf{w}}  = \\arg \\min_{\\textbf{w}} J(\\textbf X, \\textbf{y}, \\textbf{w}) =  \\arg \\min_{\\textbf{w}}\\ (C\\sum_{i=1}^{\\ell} \\log (1 + \\exp^{-y_i\\textbf{w}^{\\text{T}}\\textbf{x}_\\text{i}})+ |\\textbf{w}|^2)$$ \n",
    "\n",
    "Далее рассмотрим пример, позволяющий интуитивно понять один из смыслов регуляризации. "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
